3 Data

This analysis uses:

  • our smart meter data, downloaded via the Octopus API
  • the NG-ESO half-hourly carbon intensity data

3.1 Octopus smart meter data

See below

3.2 NG-ESO half-hourly carbon intensity data

We save this locally, if it is older than 1 day we re-download.

# this needs to be more clever - only download dates we want from API?
esoF <- here::here("data", "latest_df_fuel_ckan.csv")
if(file.exists(paste0(esoF, ".gz"))){
  message("We already have a version saved to: ", paste0(esoF, ".gz"))
  message("Loading it...")
  ngeso_dt_orig <- data.table::fread(paste0(esoF, ".gz"))
} else {
  message("We don't already have a version, downloading and saving to: ", esoF)
  ngeso_dt_orig <- data.table::fread("https://data.nationalgrideso.com/backend/dataset/88313ae5-94e4-4ddc-a790-593554d8c6b9/resource/f93d1835-75bc-43e5-84ad-12472b180a98/download/df_fuel_ckan.csv")
  # nice dateTime
  ngeso_dt_orig[, dv_start := lubridate::as_datetime(DATETIME)]
  data.table::fwrite(ngeso_dt_orig, esoF)
  dkUtils::gzipIt(esoF)
  data.table::fwrite(ngeso_dt_orig, "~/Dropbox/data/UK_NGESO/genMix/latest_df_fuel_ckan.csv") # save locally for future re-use
  dkUtils::gzipIt("~/Dropbox/data/UK_NGESO/genMix/latest_df_fuel_ckan.csv")
}
## We already have a version saved to: /Users/ben/repos/github/dataknut/octopusAPI/data/latest_df_fuel_ckan.csv.gz
## Loading it...
# if older than 12 hours, reload
maxGRID <- lubridate::as_date(max(ngeso_dt_orig$dv_start))
today <- lubridate::today()

if((today - maxGRID) > 1){
  # old data, reload
  message("But the version we have dates from ", maxGRID, " (", today - maxGRID ," days ago), downloading latest...")
  ngeso_dt_orig <- data.table::fread("https://data.nationalgrideso.com/backend/dataset/88313ae5-94e4-4ddc-a790-593554d8c6b9/resource/f93d1835-75bc-43e5-84ad-12472b180a98/download/df_fuel_ckan.csv")
  # nice dateTime
  ngeso_dt_orig[, dv_start := lubridate::as_datetime(DATETIME)]
  data.table::fwrite(ngeso_dt_orig, esoF)
  dkUtils::gzipIt(esoF)
}

ngeso_dt_orig[, dv_date := lubridate::date(dv_start)]
ngeso_dt_orig[, dv_hms := hms::as_hms(dv_start)]

# we think renewable is wind + solar, low carbon includes nuclear
#ggplot2::ggplot(t, aes(x = RENEWABLE_perc, y = LOW_CARBON_perc)) +
#  geom_point()

4 Test without authentication

Check the response code. This seems to generate errors (sometimes).

# test
url <- "https://api.octopus.energy/v1/products"
message("Getting: ", url)
## Getting: https://api.octopus.energy/v1/products
resp <- httr::GET(url)

message("Status code: ", resp$status_code)
## Status code: 200
df <- jsonlite::parse_json(resp, simplifyVector = TRUE)
## No encoding supplied: defaulting to UTF-8.
makeFlexTable(head(df$results), cap = "Example products list (first 6 rows)")

5 Tests with authentication

5.1 Basic info

url <- paste0("https://api.octopus.energy/v1/accounts/", apiParams$accountNo , "/")

resp <- httr::GET(url = url, authenticate(user = apiParams$key, password = ""))

df <- jsonlite::parse_json(resp, simplifyVector = TRUE)
## No encoding supplied: defaulting to UTF-8.
props <- data.table::as.data.table(df$properties)
makeFlexTable(head(props[, .(town, county, 
                             electricity_meter_points,gas_meter_points)]), cap = "Properties linked to this account (non-disclosive data)")

List the electricity meters by MPAN. There should be only one…

# this is a list of n mpans
length(props$electricity_meter_points)
## [1] 2
message("n MPANS listed: ", length(df$properties$electricity_meter_points))
## n MPANS listed: 2
for(n in 1:length(df$properties$electricity_meter_points)){
  print(props$electricity_meter_points[n])
}
## [[1]]
##            mpan profile_class consumption_standard
## 1 2000006198482             1                 2827
##                          meters
## 1 L78C64517, 01, STANDARD, TRUE
##                                                                                                                                                      agreements
## 1 E-1R-SUPER-GREEN-12M-20-09-22-H, E-1R-SUPER-GREEN-12M-20-09-22-H, 2020-11-01T00:00:00Z, 2020-11-21T00:00:00Z, 2020-11-01T00:00:00Z, 2021-06-30T00:00:00+01:00
##   is_export
## 1     FALSE
## 
## [[1]]
##            mpan profile_class consumption_standard
## 1 1050002522164             0                 2900
## 2 1050001805886             1                 3774
##                          meters
## 1                          NULL
## 2 19L3027004, 1, STANDARD, TRUE
##                                                                                                                                                                                                   agreements
## 1                                                                                                                                                                                                       NULL
## 2 E-1R-SUPER-GREEN-12M-20-09-22-A, E-1R-LOYAL-FIX-12M-21-10-07-A, E-1R-VAR-22-10-01-A, 2021-07-01T00:00:00+01:00, 2021-11-21T00:00:00Z, 2022-11-21T00:00:00Z, 2021-11-21T00:00:00Z, 2022-11-21T00:00:00Z, NA
##   is_export
## 1      TRUE
## 2     FALSE

List the gas meters by MPRN. There should be only one…

length(df$properties$gas_meter_points)
## [1] 2
message("n MPRNS listed: ", length(df$properties$gas_meter_points))
## n MPRNS listed: 2
for(n in 1:length(df$properties$gas_meter_points)){
  print(df$properties$gas_meter_points[n])
}
## [[1]]
##         mprn consumption_standard         meters
## 1 4256845702                15240 G4A01559730801
##                                                                                                                                                      agreements
## 1 G-1R-SUPER-GREEN-12M-20-09-22-H, G-1R-SUPER-GREEN-12M-20-09-22-H, 2020-11-01T00:00:00Z, 2020-11-21T00:00:00Z, 2020-11-01T00:00:00Z, 2021-06-30T00:00:00+01:00
## 
## [[1]]
##         mprn consumption_standard                                       meters
## 1 7825700304                15129 E6S12725512161, E6S17944211961, NOTINSTALLED
##                                                                                                                                agreements
## 1 G-1R-LOYAL-FIX-12M-21-10-07-A, G-1R-VAR-22-04-02-A, 2021-10-04T00:00:00+01:00, 2022-10-04T00:00:00+01:00, 2022-10-04T00:00:00+01:00, NA

5.2 Electricity

5.2.1 ‘Consumption’ (aka use)

See: https://www.guylipman.com/octopus/api_guide.html#s3

Start data extraction from 1st Jan 2022 as we had smart meter (re)installed in January.

Check missing dates and adjust “&page_size=100000” if required

url <- paste0("https://api.octopus.energy/v1/electricity-meter-points/", 
              apiParams$elec_import_mpan , "/",
              "meters/",
              apiParams$elec_import_serial, "/",
              "consumption/",
              "?period_from=2022-01-01T00:00Z",
              "&page_size=100000") # make sure this is large enough!
# get data via httr ----
resp <- httr::GET(url = url, authenticate(user = apiParams$key, password = ""))
df <- jsonlite::parse_json(resp, simplifyVector = TRUE) # creates a df of which 'results' = the data
## No encoding supplied: defaulting to UTF-8.
elecCons_dt <- data.table::as.data.table(df$results) # convert to dt

# derived variables ----
elecCons_dt <- makeDerivedVars(elecCons_dt)

maxTime <- max(elecCons_dt$dv_start)

hoursAgo <- now() - maxTime

# meter is SMETS2
elecCons_dt[, consumption_kWh := consumption] # for clarity - see https://developer.octopus.energy/docs/api/#list-consumption-for-a-meter

The data used here is up to 2022-12-14 23:30:00, which is 18.8 hours ago. In general the Octopus API seems to have data up to midnight last night.

5.2.2 Half-hourly analysis

Figure 5.1 shows half-hourly electricity import (‘consumption’) for the current year. Spot the power cuts…

To do: mark weekends somehow

ggplot2::ggplot(elecCons_dt, aes(x = dv_date, y = dv_hms, fill = consumption_kWh)) +
  geom_tile() +
  theme(legend.position = "bottom") +
  scale_fill_viridis_c(name = "Electricity import (kWh)") +
  labs(x = "Date",
       y = "Half-hour")
Half hourly electricity import (current year)

Figure 5.1: Half hourly electricity import (current year)

Repeat but with just the last 7 days of data - useful for checking recent appliance use and offspring effects. We do quite a lot of batch cooking on Sunday nights…

today <- lubridate::today()
plotDT <- elecCons_dt[dv_date >= today - 7]
p <- ggplot2::ggplot(plotDT, aes(x = dv_date, y = dv_hms, fill = consumption_kWh)) +
  geom_tile() +
  theme(legend.position = "bottom") +
  scale_fill_viridis_c(name = "Electricity import (kWh)") +
  labs(x = "Date",
       y = "Half-hour")

plotly::ggplotly(p)

Figure 5.2: Half hourly electricity import (current year, last 7 days) 

plotDT[, dow := lubridate::wday(dv_date, label = TRUE)]
recent_dt <- plotDT[, .(sum_kWh_elec = sum(consumption_kWh)), keyby = .(dv_date, dow, dv_peakPeriod)]
daily_totals <- plotDT[, .(Total = sum(consumption_kWh)), keyby = .(dv_date, dow)]

t <- dcast(recent_dt, dv_date + dow ~ dv_peakPeriod, val.var = sum_kWh_elec)
## Using 'sum_kWh_elec' as value column. Use 'value.var' to override
# add totals
t <- t[daily_totals]
makeFlexTable(t, digits = 2,
              cap = "Recent electricity use")
daily_totals[, fuel := "elec"]  # for future use

5.2.3 Winter 2022 SavingsSessions

Can we flex any of the above for £3/kWh?

Instigated by UK National gird, implemented (in our case) by Octopus- https://twitter.com/SavingSessions

5.2.3.1 15th Nov 2022

17:00 - 18:00 £2.50/kWh

Figure 5.3 shows how we did in the first #SavingSessions compared to the last n similar days of the week over the last 2 months.

NB: this is not the Octopus algorithm, they use the previous usage across all days (?)

startDateTime <- "2022-11-15 17:00:00" # the half-hour it starts
endDateTime <- "2022-11-15 17:30:00" # the half-hour it ends (makes data selection easier)

make_kwhComparisonPlot <- function(dt, # half hourly smart meter data
                                   startDateTime, endDateTime){
  res <- list() # results holder
  
  sessionDay <- lubridate::wday(startDateTime, label = TRUE) # you'll see why
  sessionDate <- lubridate::date(endDateTime) # you'll see why
  
  compare_HalfHourlyDT <- dt[dv_date != sessionDate & # not the session day
                               lubridate::wday(dv_date, label = TRUE) == sessionDay & # same day of the week
                               dv_date > sessionDate - 60 # 2 months
  ] # to check
  compare_HalfHourlyDT[, wday := lubridate::wday(dv_date, label = TRUE)]
  nComparisonDays <- uniqueN(compare_HalfHourlyDT$dv_date)
  
  res$nCompDays <- nComparisonDays
  
  recent_HalfHourlyDT <- compare_HalfHourlyDT[,
                                              .(mean_kWh_elec = mean(consumption_kWh)), 
                                              keyby = .(hms = dv_hms,
                                                        wday)] # to check
  recent_HalfHourlyDT[, legend_lab := paste0("Comparison mean kWh")]
  
  session_HalfHourlyDT <- elecCons_dt[dv_date == lubridate::date(startDateTime),  #  the session day
                                      .(mean_kWh_elec = mean(consumption_kWh)), 
                                      keyby = .(hms = dv_hms,
                                                wday = lubridate::wday(dv_date, label = TRUE))] # to check
  session_HalfHourlyDT[, legend_lab := "Saving session kWh"]
  plotDT <- rbind(recent_HalfHourlyDT, session_HalfHourlyDT)
  
  plotDT[, adjusted_hms := hms::as_hms(hms + (15*60))]
  
  periodAlpha <- 0.3 #  shaded rects on plots
  periodFill <- "grey50"
  ymax <- max(plotDT$mean_kWh_elec)
  ymin <- min(plotDT$mean_kWh_elec)
  xmin <- hms::as_hms(lubridate::as_datetime(startDateTime))
  xmax <- hms::as_hms(lubridate::as_datetime(endDateTime)) + 30*60 # to allow for the start time
  
  t <- plotDT[hms >= xmin & hms < xmax]
  
  wt <- dcast(t[, .(hms, wday, mean_kWh_elec, legend_lab)], hms ~ legend_lab, value.var = "mean_kWh_elec")
  
  wt[, kwh_diff := `Saving session kWh` - `Comparison mean kWh`]
  res$wt <- wt[, pc_diff := 100*(kwh_diff/`Comparison mean kWh`)]
  
  label <- paste0("SavingSession kWh (", 
                  as.Date(startDateTime),
                  ")"
                  )
  plotDT[, legend_lab := ifelse(legend_lab == "Saving session kWh",
                                label,
                                legend_lab)]
  
  res$p <- ggplot2::ggplot(plotDT, aes(x = adjusted_hms, y = mean_kWh_elec, 
                                       colour = legend_lab)) +
    geom_line() +
    geom_point() +
    annotate("rect", xmin = xmin,
             xmax = xmax,
             ymin = ymin, ymax = ymax,
             alpha = periodAlpha, fill = periodFill) +
    scale_color_discrete(name = "Legend") +
    theme(legend.position = "bottom") +
    labs(x = "Time of day",
         y = "Mean kWh per half-hour",
         caption = paste0("Data: @OctopusEnergy\nPlot: @dataknut\nPoints plotted at middle of half-hour for clarity\nNumber of number of comparison days: ", res$nCompDays)
    )
  
  return(res)
}

res <- make_kwhComparisonPlot(elecCons_dt, startDateTime, endDateTime)
res$p # the plot
First #SavingSessions - how did we do?

Figure 5.3: First #SavingSessions - how did we do? 

makeFlexTable(res$wt[, hms:= as.factor(hms)], 
              cap = paste0("kWh comparisons (hms = half hour period start, number of comparison days = ", res$nCompDays,")"),
              digits = 3) # the table

What do we conclude?

For a complicated set of pre-Xmas prep reasons (and poor institutional memory) our electricity usage in the late afternoon was much higher than usual on the session day. As a result our actions only brought usage down to ‘average’ for this time on a Tuesday (so we won’t win any points). But as a result our usage was probably way lower than it otherwise would have been judging by the spikes before and after. So if it really had been a critical peak event, we would have been helping the system but not being rewarded…

For fun - Figure 5.4 shows the mean carbon intensity in g CO2(e) per kWh for each half-hour.

make_ciComparisonPlot <- function(dt, startDateTime, endDateTime){
  res <- list() # results list
  
  sessionDay <- lubridate::wday(startDateTime, label = TRUE) # you'll see why
  sessionDate <- lubridate::date(endDateTime) # you'll see why
  
  compare_HalfHourlyCIDT <- dt[dv_start != sessionDate & # not the session day
                                 lubridate::wday(dv_date, label = TRUE) == sessionDay & # same day of the week
                                 dv_date > sessionDate - 60 # 2 months
  ] # to check
  compare_HalfHourlyCIDT[, wday := lubridate::wday(dv_date, label = TRUE)]
  nComparisonDays <- uniqueN(compare_HalfHourlyCIDT$dv_date)
  
  recent_HalfHourlyCIDT <- compare_HalfHourlyCIDT[,
                                                  .(mean_CI = mean(CARBON_INTENSITY)), 
                                                  keyby = .(hms = dv_hms,
                                                            wday)] # to check
  
  recent_HalfHourlyCIDT[, legend_lab := "Comparison mean kg CO2/kWh"]
  
  session_HalfHourlyCIDT <- dt[dv_date == lubridate::date(startDateTime),  #  the session day
                               .(mean_CI = mean(CARBON_INTENSITY)), 
                               keyby = .(hms = dv_hms,
                                         wday = lubridate::wday(dv_date, label = TRUE))] # to check
  session_HalfHourlyCIDT[, legend_lab := paste0("SavingSessions kg CO2/kWh")]
  
  plotDT <- rbind(recent_HalfHourlyCIDT, session_HalfHourlyCIDT)
  
  plotDT[, adjusted_hms := hms::as_hms(hms + (15*60))] # ad5 15 minutes so plots at middle of half-hour for clqrity
  
  periodAlpha <- 0.3 #  shaded rects on plots
  periodFill <- "grey50"
  ymax <- max(plotDT$mean_CI)
  ymin <- min(plotDT$mean_CI)
  xmin <- hms::as_hms(lubridate::as_datetime(startDateTime))
  xmax <- hms::as_hms(lubridate::as_datetime(endDateTime)) + 30*60 # to allow for the start time
  
  res$p <- ggplot2::ggplot(plotDT, aes(x = adjusted_hms, y = mean_CI, 
                                       colour = legend_lab)) +
    geom_line() +
    geom_point() +
    annotate("rect", xmin = xmin,
             xmax = xmax,
             ymin = ymin, ymax = ymax,
             alpha = periodAlpha, fill = periodFill) +
    scale_color_discrete(name = "Legend") +
    theme(legend.position = "bottom") +
    labs(x = "Time of day",
         y = "Mean carbon intensity (g CO2e/kWh)",
         caption = paste0("Data: @OctopusEnergy & @NationalGridESO \nPlot: @dataknut\nPoints plotted at middle of half-hour for clarity\nNumber of comparator days = ", nComparisonDays)
    )
  return(res)
}

res <- make_ciComparisonPlot(ngeso_dt_orig, startDateTime, endDateTime)
res$p
First #SavingSessions - comparison of carbon intensity

Figure 5.4: First #SavingSessions - comparison of carbon intensity

5.2.3.2 22nd Nov 2022

skipped this one

5.2.3.3 30th Nov 2022

17:30 - 18:30 £2.50/kWh

Figure 5.5 shows how we did in the third #SavingSessions compared to the last n similar days of the week.

startDateTime <- "2022-11-30 17:30:00" # the half-hour it starts
endDateTime <- "2022-11-30 18:00:00" # the half-hour it ends (makes data selection easier)

res <- make_kwhComparisonPlot(elecCons_dt, startDateTime, endDateTime)
res$p # the plot
Third #SavingSessions - how did we do?

Figure 5.5: Third #SavingSessions - how did we do? 

makeFlexTable(res$wt[, hms:= as.factor(hms)], 
              cap = paste0("kWh comparisons (hms = half hour period start, number of comparison days = ", res$nCompDays,")"),
              digits = 3) # the table

What do we conclude?

Cost centre #3 rushed home to get a toastie & his laundry in on the cycle before 17:00 and then went to rugby training

Cost centre #4 just went to the gym

Our dinner guest was happy to delay cooking & eating slightly

For fun - Figure 5.6 shows the mean carbon intensity in g CO2(e) per kWh for each half-hour.

res <- make_ciComparisonPlot(ngeso_dt_orig, startDateTime, endDateTime)
res$p
Third #SavingSessions - comparison of carbon intensity

Figure 5.6: Third #SavingSessions - comparison of carbon intensity

5.2.3.4 1st Dec 2022

17:00 - 18:00 £2.50/kWh

Figure 5.7 shows how we did in this #SavingSession compared to the last n similar days of the week.

startDateTime <- "2022-12-01 17:00:00" # the half-hour it starts
endDateTime <- "2022-12-01 17:30:00" # the half-hour it ends (makes data selection easier)

res <- make_kwhComparisonPlot(elecCons_dt, startDateTime, endDateTime)
res$p # the plot
Fourth #SavingSessions - how did we do?

Figure 5.7: Fourth #SavingSessions - how did we do? 

makeFlexTable(res$wt[, hms:= as.factor(hms)], 
              cap = paste0("kWh comparisons (hms = half hour period start, number of comparison days = ", res$nCompDays,")"),
              digits = 3) # the table

What do we conclude?

Cost centre #3 didn’t come home at all until 22:00

Cost centre #4 went to the gym and cooked with Cost centre #3 at 22:00

Our 2nd dinner guest of the week was also happy to delay cooking & eating slightly!

For fun - Figure 5.8 shows the mean carbon intensity in g CO2(e) per kWh for each half-hour.

res <- make_ciComparisonPlot(ngeso_dt_orig, startDateTime, endDateTime)
res$p
Fourth #SavingSessions - comparison of carbon intensity

Figure 5.8: Fourth #SavingSessions - comparison of carbon intensity

5.2.3.5 12th Dec 2022

17:00 - 18:00 £2.50/kWh

Figure 5.9 shows how we did in this #SavingSession compared to the last n similar days of the week.

startDateTime <- "2022-12-12 17:00:00" # the half-hour it starts
endDateTime <- "2022-12-12 18:30:00" # the half-hour it ends (makes data selection easier)

res <- make_kwhComparisonPlot(elecCons_dt, startDateTime, endDateTime)
res$p # the plot
Fifth #SavingSessions - how did we do?

Figure 5.9: Fifth #SavingSessions - how did we do? 

makeFlexTable(res$wt[, hms:= as.factor(hms)], 
              cap = paste0("kWh comparisons (hms = half hour period start, number of comparison days = ", res$nCompDays,")"),
              digits = 3) # the table

What do we conclude?

Pre-planned trip to Grandma’s worked well…

For fun - Figure 5.10 shows the mean carbon intensity in g CO2(e) per kWh for each half-hour.

res <- make_ciComparisonPlot(ngeso_dt_orig, startDateTime, endDateTime)
res$p
Fifth #SavingSessions - comparison of carbon intensity

Figure 5.10: Fifth #SavingSessions - comparison of carbon intensity

5.2.4 Daily analysis

Figure 5.11 shows the total daily kWh import with a smoothed curve for each day as shown.

plotDT <- elecCons_dt[, .(sum_kWh = sum(consumption_kWh),
                          mean_kWh = mean(consumption_kWh),
                          nObs = .N), keyby = .(dv_date, dv_weekend)]

ggplot2::ggplot(plotDT, aes(x = dv_date, y = sum_kWh, 
                            colour = dv_weekend)) +
  geom_point() +
  geom_smooth() +
  #facet_grid(dv_peakPeriod ~ .) +
  scale_colour_viridis_d(name = "Weekend") +
  labs(x = "Date",
       y = "Daily kWh")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Mean half-hourly electricity import per day (kWh, current year)

Figure 5.11: Mean half-hourly electricity import per day (kWh, current year)

Figure 5.12 shows the mean daily kWh import with a smoothed curve for each period as defined below. The periods do not have the same number of half-hours so we use the mean as a comparator.

Early morning is effectively our baseload.

# check periods
t <- elecCons_dt[, .(min = min(dv_hms),
                     max = max(dv_hms)),
                 keyby = .(dv_peakPeriod)]
t
##    dv_peakPeriod      min      max
## 1: Early morning 00:00:00 06:30:00
## 2:  Morning peak 07:00:00 08:30:00
## 3:      Day time 09:00:00 15:30:00
## 4:  Evening peak 16:00:00 19:30:00
## 5:  Late evening 20:00:00 23:30:00
t <- elecCons_dt[, .(mean_kWh = mean(consumption_kWh),
                     sd_kWh = sd(consumption_kWh),
                     min_kWh = min(consumption_kWh),
                     max_kWh = max(consumption_kWh)),
                 keyby = .(dv_peakPeriod)]
makeFlexTable(t, digits = 3, cap = "Summary stats (all half-hourly data)")
plotDT <- elecCons_dt[, .(sum_kWh = sum(consumption_kWh),
                          mean_kWh = mean(consumption_kWh),
                          nObs = .N), keyby = .(dv_date, dv_peakPeriod)]
totalDT <- elecCons_dt[, .(sum_kWh = sum(consumption_kWh),
                           mean_kWh = mean(consumption_kWh),
                           nObs = .N), keyby = .(dv_date)]
totalDT[, dv_peakPeriod := "All periods"]

plotDT <- rbind(plotDT, totalDT)

ggplot2::ggplot(plotDT, aes(x = dv_date, y = mean_kWh, 
                            colour = dv_peakPeriod)) +
  geom_line() +
  geom_smooth() +
  #facet_grid(dv_peakPeriod ~ .) +
  theme(legend.position = "bottom") +
  guides(colour = guide_legend (ncol = 3)) +
  scale_colour_viridis_d(name = "Peak period") +
  labs(x = "Date",
       y = "Mean kWh per half-hour")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Mean half-hourly electricity import (kWh) by peak period (current year)

Figure 5.12: Mean half-hourly electricity import (kWh) by peak period (current year)

What’s happening here?

  • we have PV (but not a lot) -> summer Day time
  • we switched from a gas hob to induction in July -> evening peak
  • we got rid of a freezer with a bad seal in July (-> early morning baseload) and also swapped a built-in fridge-freezer for a stand-alone fridge (& freezer)
  • we set the hot water circulation loop to ‘just in time’ (morning peak). Well, kind of. Might affect the gas more.

5.2.5 PV generation

Not Octopus data - irregular reads by us. Load from .xslsx

Figure 5.13 compares our PV generation with overall grid import. These are mean values. We are occasionally exporting (Table ?? has zeros) but we do not (currently) have access to the export data.

library(openxlsx)
f <- "~/Dropbox/Home/houses/whiteHorseMews/2WhiteHorseMewsRunningCosts.xlsx"
pvGen <- openxlsx::read.xlsx(f,
                             sheet = "PV",
                             detectDates = TRUE)
pvGen_DT <- data.table::as.data.table(pvGen)
pvGen_DT[, meankWhGen := diff/n.days]
pvGen_DT[, dv_date := lubridate::as_datetime(Date)]
pvGen_DT[, month := lubridate::month(dv_date)]

monthlyPvGen_DT <- pvGen_DT[, .(hh_meankWh = mean(meankWhGen)/48 # input value is mean per day
),
keyby = .(month = lubridate::month(dv_date),
          year = lubridate::year(dv_date))
]
monthlyPvGen_DT[, value := "PV generation"]

monthlyElec_DT <- elecCons_dt[, .(hh_meankWh = mean(consumption_kWh)),
                              keyby = .(month = lubridate::month(dv_date),
                                        year = lubridate::year(dv_date))]
monthlyElec_DT[, value := "Grid import"]

plotDT <- rbind(monthlyPvGen_DT, monthlyElec_DT)

ggplot2::ggplot(plotDT, aes(x = lubridate::month(month, label = TRUE), y = hh_meankWh, 
                            colour = value, group = value)) +
  geom_line()+
  scale_color_discrete(name = "Type")+
  facet_grid(year ~ .) +
  labs(x = "Month")
## Warning: Removed 1 row containing missing values (`geom_line()`).
Compare grid import and PV gen

Figure 5.13: Compare grid import and PV gen

To do: model value of PV gen if we were to use all of it.

5.2.6 PV Export

This will be a new MPAN but specified as export - although the url will still say consumption. We do not have this even though the PV is exporting on (some) days.

It may be that we only get this data if we sign up for the export tariff.

See https://www.guylipman.com/octopus/api_guide.html#s3

5.2.7 Electricity emissions

In theory our emissions from electricity use are zero because we are on a renewable-only tariff. But life is not so simple. We don’t have a private wire to a wind turbine so the electrons we import (stick with it) are as averagely green as all the rest.

We also con’t avoid the ‘Well To Tank’ emissions and those associated with transmission losses.

To further complicate things there are at least two different ways to estimate our emissions.

  • use the annual BEIS emissions factor and multiply by the kWh in question - be it half-hourly, daily, annual, whatever. That’s the simple way.
  • use the NG-ESO half-hourly emissions factors which reflect the generation mix (with some caveats) of the grid at half-hourly intervals

Does it matter? you cry. Well it might. If we’ve been able to ‘flex’ our usage in line with @theBakingForecast then who knows, maybe we’ll be concentrating our usage in times when the grid is actually drawing on more renewables.

So let’s take a look. We’ll do both the BEIS-based and NG-ESO based calculations to see. For now we’ll ignore the WTT and the T&D losses to keep the results comparable. We’ll come back to that later.

rmdParams$BEIS_elec_ci <-   0.21233

For the BEIS method, we’ll have to use the 2021 emissions factor as the 2022 value is not yet available.

For 2021 this is: 0.21233 Kg CO2e/kWh

For the NG-ESO method we use the NG-ESO half-hourly carbon intensity data that match to the half-hours in our electricity use dataset.

Mean half-hourly carbon intensity from the NG-ESO data for the data period was NA Kg CO2e/kWh. If this is substantially different to the BEIS 2021 value of 0.2123 Kg CO2e/kWh, we would expect emissions estimates using the NG-ESO factor to differ.

# merge to usage data
setkey(ngeso_dt_orig, dv_start)
setkey(elecCons_dt, dv_start)
elecCons_dt <- ngeso_dt_orig[, .(dv_start, CARBON_INTENSITY, LOW_CARBON_perc, RENEWABLE_perc)][elecCons_dt] # keeps match to our electricity use

# there will be NAs if some datetimes are missing from ngeso_dt_orig

For context, Figure 5.14 summarises the mean half-hourly carbon intensity by month for the data period. We can clearly see that February 2022 was a very low carbon month… in fact it was a very windy month with 3 named storms.

elecCons_dt[, dv_month := lubridate::month(dv_date, label = TRUE)]
ggplot2::ggplot(elecCons_dt, aes(x = dv_month, y = CARBON_INTENSITY)) +
  geom_violin(draw_quantiles = c(0.5)) +
  #geom_jitter() +
  #geom_boxplot() +
  labs(x = "Month",
       y = "Half-hourly carbon intensity",
       caption = "Median drawn")
## Warning: Removed 23 rows containing non-finite values (`stat_ydensity()`).
Monthly mean carbon intensity for the data period by month (NG-ESO data)

Figure 5.14: Monthly mean carbon intensity for the data period by month (NG-ESO data)

Figure 5.15 shows our half-hourly electricity kWh use vs halfhourly carbon intensity. Ideally we want a negative correlation showing that we use the most electricity when it is ‘greenest’ (carbon intensity is lowest). Doesn’t look too good, aye?

ggplot2::ggplot(elecCons_dt, aes(x = CARBON_INTENSITY, y = consumption_kWh, colour = RENEWABLE_perc)) +
  geom_point() +
  facet_wrap(. ~ dv_peakPeriod) +
  geom_smooth() +
  scale_color_continuous(name = "% renewables", low = "red", high = "green") +
  theme(legend.position = "bottom") +
  labs(y = "Halfhourly electricity kWh")
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
## Warning: Removed 23 rows containing non-finite values (`stat_smooth()`).
## Warning: The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
## Warning: Removed 23 rows containing missing values (`geom_point()`).
Own half-hourly electricity kWh vs NG-ESO halfhourly carbon intensity

Figure 5.15: Own half-hourly electricity kWh vs NG-ESO halfhourly carbon intensity

What if we visualise using a box plot according to carbon intensity decile? So this means we divide the carbon intensity values into 10 equal groups - deciles. This is Figure 5.16. Doesn’t look too good either - our median usage (the horizontal bar in the boxes) seems to trend slightly upwards as we move to higher carbon intensity deciles.

# this will generate NAs if the CI data is missing for some of the (most recent) dateTimes 
elecCons_dt[, CI_deciles := cut_number(CARBON_INTENSITY, n = 10)]

ggplot2::ggplot(elecCons_dt, aes(x = CI_deciles, y = consumption_kWh)) +
  geom_boxplot() +
  labs(y = "Halfhourly electricity kWh",
       x = "Carbon intensity decile")
Boxing clever

Figure 5.16: Boxing clever

So what if we just add up all our electricity kWh by carbon intensity decile? Do we use more low carbon kWh? This is Figure 5.17. Nah. The bakingforecast isn’t going to like us…

t <- elecCons_dt[, .(sumkWh = sum(consumption_kWh, na.rm = TRUE),
                     meankWh = mean(consumption_kWh, na.rm = TRUE)),
                 keyby = .(CI_deciles)]

ggplot2::ggplot(t, aes(x = CI_deciles, y = sumkWh)) +
  geom_col() +
  labs(y = "Sum kWh",
       x = "Carbon intensity decile")
Sum of electricity kWh by carbon intensity decile

Figure 5.17: Sum of electricity kWh by carbon intensity decile

Out of interest, do our emissions values look very different if we apply the BEIS 2021 annual factor to our total electricity kWh to date compared to applying the NG-ESO half-hourly values?

elecCons_dt[, KgCO2_ngeso := consumption_kWh * (CARBON_INTENSITY/1000)] # convert to kg
t <- elecCons_dt[, .(sumkWh = sum(consumption_kWh),
                     sumKgCO2_ngeso = sum(KgCO2_ngeso, na.rm = TRUE))]

t[, sumKgCO2_beis :=  sumkWh * rmdParams$BEIS_elec_ci]

makeFlexTable(t, cap = "Comparing emissions estimation methods using electricity kWh to date")
t <- elecCons_dt[, .(sumkWh = sum(consumption_kWh),
                     sumKgCO2_ngeso = sum(KgCO2_ngeso, na.rm = TRUE)),
                 keyby = .(dv_month)]

t[, sumKgCO2_beis :=  sumkWh * rmdParams$BEIS_elec_ci]

plotDT <- melt(t, id.vars = "dv_month")

ggplot2::ggplot(plotDT[ variable != "sumkWh",], aes(x = dv_month, y = value, fill = variable)) +
  geom_col(position = "dodge") +
  scale_color_discrete(name = "Method") +
  labs(y = "Kg CO2",
       x = "Month")

As we’d expect from the comparison of the values above, Table 5.8 suggests that it does. In fact our ‘in use’ NG-ESO based emissions are 88.48, 63.96, 91.43, 87.07, 90.3, 88.01, 99.7, 106.04, 100.82, 76.72, 83.61, 115.14 % of our BEIS-based emissions depending on the month in question.

If we compare the monthly values we can see the biggest difference was in February, a month we have already identified as being more ‘low carbon’ (see Figure 5.14).

5.2.8 Electricity costs (estimates)

Analyse costs using:

The latter are slightly different from the assumed to be at UK price cap: £0.34 / kWh & £0.46 ( see Ofgem)

prices <- readxl::read_xlsx(here::here("data", "prices.xlsx"))
pricesDT <- data.table::as.data.table(prices)

Yes, I know I can extract our exact tariff from the octopus API…

daily_elec <- elecCons_dt[, .(sum_kWh = sum(consumption_kWh),
                           nObs = .N), keyby = .(dv_date)]

# extract from pricesDT
# must be an easier way

daily_elec[, kwh_p := ifelse(dv_date < lubridate::as_date("2022-11-21"),
                               pricesDT[fuel == "elec" & component == "kWh" &
                  dateEnd == "2022-11-21", price],
                               pricesDT[fuel == "elec" & component == "kWh" &
                  dateStart == "2022-11-22", price])]
daily_elec[, sc_p := ifelse(dv_date < lubridate::as_date("2022-11-21"),
                               pricesDT[fuel == "elec" & component == "sc" &
                  dateEnd == "2022-11-21", price],
                               pricesDT[fuel == "elec" & component == "sc" &
                  dateStart == "2022-11-22", price])]

daily_elec[, cost := ((sum_kWh * kwh_p) + sc_p)]
daily_elec[, month := lubridate::month(dv_date, label = TRUE)]

ggplot2::ggplot(daily_elec, aes(x = dv_date, y = cost, colour = month)) +
  geom_line() +
  geom_smooth() +
  geom_vline(xintercept = lubridate::as_date("2022-11-21")) +
  labs(y = "Daily cost £",
       caption = "Tariff change/price cap/EPG 2022 date shown\nSmoothed within month")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Daily electricity costs

Figure 5.18: Daily electricity costs 

daily_elec[, month_floor := lubridate::floor_date(dv_date, "months")]
monthly_elec <- daily_elec[, .(sum_kWh = sum(sum_kWh),
                               cost = sum(cost)),
                           keyby = .(month_floor)]

ggplot2::ggplot(monthly_elec, aes(x = month_floor, y = cost)) +
  geom_col() +
  geom_vline(xintercept = lubridate::as_date("2022-11-21")) +
  labs(y = "Monthly cost £",
       x = "Month", 
       caption = "Tariff change/price cap/EPG 2022 date shown\nBeware incomplete months")
Monthly electricity costs

Figure 5.19: Monthly electricity costs 

message("Projected annual elec total kWh")
## Projected annual elec total kWh
projAnnual_elec_kWh <- mean(daily_elec$sum_kWh)*365
projAnnual_elec_kWh
## [1] 3739.295
message("########")
## ########
message("# Prices to 21st Nov 2022")
## # Prices to 21st Nov 2022
# get_price <- function(dt, fuel, component, dateEnd){
#   p <- dt[fuel == fuel &
#             component == component &
#             dateEnd  == dateEnd,
#           price]
#   return(p)
# }

kWh_p <- pricesDT[fuel == "elec" & component == "kWh" &
                  dateEnd == "2022-11-21", price]
sc_p <- pricesDT[fuel == "elec" & component == "sc" &
                  dateEnd == "2022-11-21", price]

message("Projected annual elec cost @ £ ", kWh_p,
        " per kWh & standing charge at £ ", sc_p,
        " per day")
## Projected annual elec cost @ £ 0.2408 per kWh & standing charge at £ 0.2401 per day
projAannualCost <- (projAnnual_elec_kWh*kWh_p)+(365*sc_p) # standing charge
projAannualCost
## [1] 988.0587
message("Mean projected monthly cost: £")
## Mean projected monthly cost: £
projAannualCost/12
## [1] 82.33823
message("########")
## ########
message("# From to 22nd Nov 2022 - prioce cap")
## # From to 22nd Nov 2022 - prioce cap
kWh_p <- pricesDT[fuel == "elec" & component == "kWh" &
                  dateStart == "2022-11-22", price]
sc_p <- pricesDT[fuel == "elec" & component == "sc" &
                  dateStart == "2022-11-22", price]

message("Projected annual elec cost @ £ ", kWh_p,
        " per kWh & standing charge at £ ", sc_p,
        " per day")
## Projected annual elec cost @ £ 0.3506 per kWh & standing charge at £ 0.3729 per day
annualCost_elecCapped <- (projAnnual_elec_kWh*kWh_p)+(365*sc_p) # standing charge
annualCost_elecCapped
## [1] 1447.105
message("Projected mean monthly £ under price cap")
## Projected mean monthly £ under price cap
monthlyCost_capped <- annualCost_elecCapped/12
monthlyCost_capped
## [1] 120.5921
message("That's an increase of ", round(100*((annualCost_elecCapped-projAannualCost)/projAannualCost),2),  " % points")
## That's an increase of 46.46 % points

TO DO: fix £ analysis to use tariff from API

URL will be something like

https://api.octopus.energy/v1/products/LOYAL-FIX-12M-21-10-07/electricity-tariffs/E-1R-LOYAL-FIX-12M-21-10-07-A/standard-unit-rates/

5.3 Gas ‘Consumption’ (aka use)

We need to convert the gas consumption from m3 to kWh - see https://developer.octopus.energy/docs/api/#list-consumption-for-a-meter

gasM3TokWh <- 11.36

We use a multiplier of 11.36 kWh/m3 (https://www.theenergyshop.com/guides/how-to-convert-gas-units-to-kwh)

Check for missing dates and adjust “&page_size=100000” if needed

url <- paste0("https://api.octopus.energy/v1/gas-meter-points/", 
              apiParams$gas_mpan , "/",
              "meters/",
              apiParams$gas_serial, "/",
              "consumption",
              "?period_from=2022-01-01T00:00Z",
              "&page_size=100000") # make sure is large enough
resp <- httr::GET(url = url, authenticate(user = apiParams$key, password = ""))
df <- jsonlite::parse_json(resp, simplifyVector = TRUE)
## No encoding supplied: defaulting to UTF-8.
gasCons_dt <- data.table::as.data.table(df$results)
gasCons_dt <- makeDerivedVars(gasCons_dt)

# gas 'consumption' is m3 - https://developer.octopus.energy/docs/api/#list-consumption-for-a-meter
# convert to kWh
gasCons_dt[, consumption_m3 := consumption]
gasCons_dt[, consumption_kWh := consumption * gasM3TokWh]

Note that this data starts later as we finally got the original un-registered smart meter replaced in February.

5.3.1 Half-hourly anlaysis

Figure 5.20 shows half-hourly gas import (‘consumption’) for the current year. The power cuts are even easier to see here. Interestingly the pattern after the gas boiler was serviced in June is more varied. What did he change?

ggplot2::ggplot(gasCons_dt, aes(x = dv_date, y = dv_hms, fill = consumption_kWh)) +
  geom_tile() +
  theme(legend.position = "bottom") +
  scale_fill_viridis_c(name = "Gas import (kWh)") +
  labs(x = "Date",
       y = "Half-hour")
Half-hourly gas consumption (current year)

Figure 5.20: Half-hourly gas consumption (current year)

Repeat but with just the last 7 days of data - useful for checking recent appliance use and offspring effects.

today <- lubridate::today()
plotDT <- gasCons_dt[dv_date >= today - 7]
p <- ggplot2::ggplot(plotDT, aes(x = dv_date, y = dv_hms, fill = consumption_kWh)) +
  geom_tile() +
  theme(legend.position = "bottom") +
  scale_fill_viridis_c(name = "Gas import (kWh)") +
  labs(x = "Date",
       y = "Half-hour")

plotly::ggplotly(p)

Figure 5.21: Half hourly gas import (current year, last 7 days) 

plotDT[, dow := lubridate::wday(dv_date, label = TRUE)]
recent_dt <- plotDT[, .(sum_kWh_elec = sum(consumption_kWh)), keyby = .(dv_date, dow, dv_peakPeriod)]
daily_totals <- plotDT[, .(Total = sum(consumption_kWh)), keyby = .(dv_date, dow)]

t <- dcast(recent_dt, dv_date + dow ~ dv_peakPeriod, val.var = sum_kWh_elec)
## Using 'sum_kWh_elec' as value column. Use 'value.var' to override
# add totals
t <- t[daily_totals]
makeFlexTable(t, digits = 2,
              cap = "Recent gas use")

5.3.2 Daily analysis

Figure ?? shows the mean daily kWh import with a smoothed curve.

To do: mark weekends etc

plotDT <- gasCons_dt[, .(sum_kWh = sum(consumption_kWh),
                         mean_kWh = mean(consumption_kWh),
                         nObs = .N), keyby = .(dv_date)]

ggplot2::ggplot(plotDT, aes(x = dv_date, y = sum_kWh)) +
  geom_point() +
  geom_smooth() +
  theme(legend.position = "bottom") +
  labs(x = "Date",
       y = "Sum kWh per day")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Daily gas consumption (current year)

Figure 5.22: Daily gas consumption (current year) 

plotDT <- gasCons_dt[, .(sum_kWh = sum(consumption_kWh),
                         mean_kWh = mean(consumption_kWh),
                         nObs = .N), keyby = .(dv_date, dv_weekend)]

ggplot2::ggplot(plotDT, aes(x = dv_date, y = sum_kWh, 
                            colour = dv_weekend)) +
  geom_point() +
  geom_smooth() +
  theme(legend.position = "bottom") +
  guides(colour = guide_legend (ncol = 3)) +
  scale_colour_viridis_d(name = "Weekend") +
  labs(x = "Date",
       y = "Sum kWh per day")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Daily gas consumption (current year)

Figure 5.23: Daily gas consumption (current year)

Repeat for mean

ggplot2::ggplot(plotDT, aes(x = dv_date, y = mean_kWh, 
                            colour = dv_weekend)) +
  geom_point() +
  geom_smooth() +
  theme(legend.position = "bottom") +
  guides(colour = guide_legend (ncol = 3)) +
  scale_colour_viridis_d(name = "Weekend") +
  labs(x = "Date",
       y = "Mean kWh per day")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Figure 5.24 shows the mean daily kWh import with a smoothed curve by period of the day.

To do: mark weekends etc

plotDT <- gasCons_dt[, .(sum_kWh = sum(consumption_kWh),
                         mean_kWh = mean(consumption_kWh),
                         nObs = .N), keyby = .(dv_date, dv_peakPeriod, dv_weekend)]

ggplot2::ggplot(plotDT, aes(x = dv_date, y = mean_kWh, 
                            colour = dv_peakPeriod)) +
  geom_line() +
  geom_smooth() +
  #facet_grid(dv_peakPeriod ~ .) +
  theme(legend.position = "bottom") +
  guides(colour = guide_legend (ncol = 3)) +
  scale_colour_viridis_d(name = "Peak period") +
  labs(x = "Date",
       y = "Mean kWh per period")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Daily gas consumption by peak period (current year)

Figure 5.24: Daily gas consumption by peak period (current year)

5.3.3 Gas emissions

This is much more simple. We can only apply the BEIS 2021 value as there are no time-varying emissions factors for gas.

rmdParams$BEIS_gas_ci <- 0.20297

As before, for the BEIS method we’ll have to use the 2021 emissions factor as the 2022 value is not yet available.

For 2021 this is: 0.20297 Kg CO2e/kWh

gasCons_dt[, KgCO2_beis := consumption_kWh * rmdParams$BEIS_gas_ci] 
t <- gasCons_dt[, .(sumkWh = sum(consumption_kWh),
                    sumKgCO2_beis = sum(KgCO2_beis))]

makeFlexTable(t, cap = "Emissions estimation using gas kWh to date")

5.3.4 Gas costs

Analyse costs using:

The latter are (currently) the same as the at UK price cap: £0.1031 / kWh & £0.2684 ( see Ofgem)

Yes, I know I can extract our exact tariff from the octopus API…

daily_gas <- gasCons_dt[, .(sum_kWh = sum(consumption_kWh),
                           nObs = .N), keyby = .(dv_date)]

# extract from pricesDT
# must be an easier way

daily_gas[, kwh_p := ifelse(dv_date < lubridate::as_date("2022-10-04"),
                               pricesDT[fuel == "gas" & component == "kWh" &
                  dateEnd %like% "2022-10-03", price], # why does it need to be like??
                               pricesDT[fuel == "gas" & component == "kWh" &
                  dateStart %like% "2022-10-04", price])]
daily_gas[, sc_p := ifelse(dv_date < lubridate::as_date("2022-10-04"),
                               pricesDT[fuel == "gas" & component == "sc" &
                  dateEnd %like% "2022-10-03", price],
                               pricesDT[fuel == "gas" & component == "sc" &
                  dateStart %like% "2022-10-04", price])]

daily_gas[, cost := ((sum_kWh * kwh_p) + sc_p)]
daily_gas[, month := lubridate::month(dv_date, label = TRUE)]

ggplot2::ggplot(daily_gas, aes(x = dv_date, y = cost, colour = month)) +
  geom_line() +
  geom_smooth() +
  geom_vline(xintercept = lubridate::as_date("2022-10-04")) +
  labs(y = "Daily cost £",
       caption = "Tariff change/price cap/EPG 2022 date shown\nSmoothed within month")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Daily gas costs

Figure 5.25: Daily gas costs 

lastWeek <- max(daily_gas$dv_date) - 7

makeFlexTable(daily_gas[dv_date > lastWeek, .(dv_date, 
                                              day = lubridate::wday(dv_date, label = TRUE),
                                              sum_kWh, nObs, cost)], digits = 2,
              cap = "Recent daily gas cost")
daily_gas[, month_floor := lubridate::floor_date(dv_date, "months")]
monthly_gas <- daily_gas[, .(sum_kWh = sum(sum_kWh),
                               cost = sum(cost)),
                           keyby = .(month_floor)]

ggplot2::ggplot(monthly_gas, aes(x = month_floor, y = cost)) +
  geom_col() +
  geom_vline(xintercept = lubridate::as_date("2022-10-04")) +
  labs(y = "Monthly cost £",
       x = "Month", 
       caption = "Tariff change/price cap/EPG 2022 date shown\nBeware incomplete months")
Monthly gas costs

Figure 5.26: Monthly gas costs 

message("Projected annual gas total kWh")
## Projected annual gas total kWh
projAnnual_gas_kWh <- mean(daily_gas$sum_kWh)*365
projAnnual_gas_kWh
## [1] 16494.85
message("########")
## ########
message("# Prices to 10th October 2022")
## # Prices to 10th October 2022
# get_price <- function(dt, fuel, component, dateEnd){
#   p <- dt[fuel == fuel &
#             component == component &
#             dateEnd  == dateEnd,
#           price]
#   return(p)
# }

kWh_p <- pricesDT[fuel == "gas" & component == "kWh" &
                  dateEnd %like% "2022-10-03", price]
sc_p <- pricesDT[fuel == "gas" & component == "sc" &
                  dateEnd %like% "2022-10-03", price]

message("Projected annual gas cost @ £ ", kWh_p,
        " per kWh & standing charge at £ ", sc_p,
        " per day")
## Projected annual gas cost @ £ 0.0719 per kWh & standing charge at £ 0.261 per day
projAannualCost <- (projAnnual_gas_kWh*kWh_p)+(365*sc_p) # standing charge
projAannualCost
## [1] 1281.245
message("Mean projected monthly cost: £")
## Mean projected monthly cost: £
projAannualCost/12
## [1] 106.7704
message("########")
## ########
message("# From 4th Oct 2022 - price cap")
## # From 4th Oct 2022 - price cap
kWh_p <- pricesDT[fuel == "gas" & component == "kWh" &
                  dateStart %like% "2022-10-04", price]
sc_p <- pricesDT[fuel == "gas" & component == "sc" &
                  dateStart %like% "2022-10-04", price]

message("Projected annual gas cost @ £ ", kWh_p,
        " per kWh & standing charge at £ ", sc_p,
        " per day")
## Projected annual gas cost @ £ 0.1031 per kWh & standing charge at £ 0.2684 per day
annualCost_gasCapped <- (projAnnual_gas_kWh*kWh_p)+(365*sc_p) # standing charge
annualCost_gasCapped
## [1] 1798.585
message("Projected mean monthly £ under price cap")
## Projected mean monthly £ under price cap
monthlyCost_capped <- annualCost_gasCapped/12
monthlyCost_capped
## [1] 149.8821
message("That's an increase of ", round(100*((annualCost_gasCapped-projAannualCost)/projAannualCost),2),  " % points")
## That's an increase of 40.38 % points

6 Total energy & costs

Monthly estimated costs (averaged over 12 months to guide monthly payment value).

These use the latest Octopus flexible tariff from Oct 2022 which are protected by the government’s Energy Price Guarantee:

t <- pricesDT[dateStart > as.Date("2022-10-01")]
t[, dateStart := as.Date(dateStart)]
makeFlexTable(t, cap = "Price cap prices", digits = 3)

According to our EPC:

  • floor area: 163 m2
  • Air permeability 5.0 m³/h.m² (as tested)
  • The primary energy use for this property per year is 65 kilowatt hours per square metre (kWh/m2) or 1.0595^{4} kWh/year (“Primary energy use is a measure of the energy required for lighting, heating and hot water in a property”)
  • This property produces 1.9 tonnes of CO2
  • Space heating 6511 kWh per year
  • Water heating 2351 kWh per year
  • total of 8862 kWh/year

First we’ll compare with hot water…

message("Assume our August-September gas use is just for hot water (almost certainly always true - no gas hob and no heating on)")
## Assume our August-September gas use is just for hot water (almost certainly always true - no gas hob and no heating on)
hw <- gasCons_dt[dv_date >= as.Date("2022-08-01") &
                                     dv_date < as.Date("2022-10-01"), .(daily_kWh = sum(consumption_kWh),
                                                                        nObs = .N),
                 keyby = .(dv_date)]

message("Mean daily gas use in this period = ", round(mean(hw$daily_kWh),2), " kWh")
## Mean daily gas use in this period = 31.92 kWh
annualHW <- 365 * mean(hw$daily_kWh)
message("So estimated annual gas for hot water = ", round(annualHW), " kWh (assuming constant hot water use all year round)")
## So estimated annual gas for hot water = 11652 kWh (assuming constant hot water use all year round)
message("That's ", round(annualHW/2351,2), " times the EPC estimate of 2351 kWh... ")
## That's 4.96 times the EPC estimate of 2351 kWh...

And now heating…

message("Assume the estimate of hot water kWh is true")
## Assume the estimate of hot water kWh is true
message("So estimated annual gas for heating = ", round(projAnnual_gas_kWh - annualHW,2), " kWh")
## So estimated annual gas for heating = 4842.68 kWh
message("That's ", round((projAnnual_gas_kWh - annualHW)/6511,2), " times the EPC estimate of 6511 kWh... ")
## That's 0.74 times the EPC estimate of 6511 kWh...

Now switching to overall costs…

message("# Total projected energy")
## # Total projected energy
message(round(projAnnual_gas_kWh + projAnnual_elec_kWh), " kWh")
## 20234 kWh
message("This is ", round(100*((projAnnual_gas_kWh + projAnnual_elec_kWh)/(65*163))), " % of our EPC 'primary energy'")
## This is 191 % of our EPC 'primary energy'
message("#####")
## #####
message("# Total costs")
## # Total costs
message("Monthly total £:", round((annualCost_gasCapped + annualCost_elecCapped)/12))
## Monthly total £:270
annualTotalcapped <- (annualCost_gasCapped + annualCost_elecCapped)

message("Annual total £:", round(annualTotalcapped,2))
## Annual total £:3245.69

If we then subtract the Energy Bill Support Scheme £400 this gives

# subtract £400 Energy Bill Support Scheme
annualTotalcappedAdjusted <- annualTotalcapped - 400
message("Annual adjusted total £:", round(annualTotalcappedAdjusted,2))
## Annual adjusted total £:2845.69
message("Monthly adjusted total £:", round(annualTotalcappedAdjusted/12,2))
## Monthly adjusted total £:237.14

Note that octopus appear to be ‘paying’ this in instalments of £67 by reducing the direct debit value.

NB - this is a simple average across all days/months and takes no account of usage trends (see above). The octopus estimator is smarter than that :-)

7 Annexes

7.1 Data descriptions

7.1.1 Electricity consumption

Use skmir::skim()to summarise.

skimr::skim(elecCons_dt)
Table 7.1: Data summary
Name elecCons_dt
Number of rows 16704
Number of columns 15
Key dv_start
_______________________
Column type frequency:
character 3
Date 1
difftime 1
factor 3
numeric 6
POSIXct 1
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
interval_start 0 1 20 25 0 16704 0
interval_end 0 1 20 25 0 16704 0
dv_weekend 0 1 6 8 0 3 0

Variable type: Date

skim_variable n_missing complete_rate min max median n_unique
dv_date 0 1 2022-01-01 2022-12-14 2022-06-23 348

Variable type: difftime

skim_variable n_missing complete_rate min max median n_unique
dv_hms 0 1 0 secs 84600 secs 42300 secs 48

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
dv_peakPeriod 0 1 FALSE 5 Ear: 4872, Day: 4872, Eve: 2784, Lat: 2784
dv_month 0 1 TRUE 12 Jan: 1488, Mar: 1488, May: 1488, Jul: 1488
CI_deciles 23 1 FALSE 10 (19: 1726, (88: 1698, (24: 1697, (15: 1694

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
CARBON_INTENSITY 23 1 185.07 64.04 39.0 136.00 193.00 237.00 323.00 ▃▅▇▇▂
LOW_CARBON_perc 23 1 52.70 14.86 18.9 40.80 51.10 63.60 91.50 ▂▇▇▅▂
RENEWABLE_perc 23 1 31.63 15.48 2.2 19.20 30.30 43.80 71.10 ▅▇▇▅▂
consumption 0 1 0.21 0.21 0.0 0.09 0.15 0.26 1.97 ▇▁▁▁▁
consumption_kWh 0 1 0.21 0.21 0.0 0.09 0.15 0.26 1.97 ▇▁▁▁▁
KgCO2_ngeso 23 1 0.04 0.05 0.0 0.01 0.03 0.05 0.53 ▇▁▁▁▁

Variable type: POSIXct

skim_variable n_missing complete_rate min max median n_unique
dv_start 0 1 2022-01-01 2022-12-14 23:30:00 2022-06-23 23:45:00 16704

7.1.2 Gas consumption

Use skmir::skim()to summarise.

skimr::skim(gasCons_dt)
Table 7.2: Data summary
Name gasCons_dt
Number of rows 14277
Number of columns 11
Key NULL
_______________________
Column type frequency:
character 3
Date 1
difftime 1
factor 1
numeric 4
POSIXct 1
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
interval_start 0 1 20 25 0 14277 0
interval_end 0 1 20 25 0 14277 0
dv_weekend 0 1 6 8 0 3 0

Variable type: Date

skim_variable n_missing complete_rate min max median n_unique
dv_date 0 1 2022-02-11 2022-12-14 2022-07-11 301

Variable type: difftime

skim_variable n_missing complete_rate min max median n_unique
dv_hms 0 1 0 secs 84600 secs 12:00:00 48

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
dv_peakPeriod 0 1 FALSE 5 Day: 4166, Ear: 4158, Eve: 2384, Lat: 2381

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
consumption 0 1 0.08 0.16 0 0 0 0.11 1.19 ▇▁▁▁▁
consumption_m3 0 1 0.08 0.16 0 0 0 0.11 1.19 ▇▁▁▁▁
consumption_kWh 0 1 0.95 1.79 0 0 0 1.26 13.53 ▇▁▁▁▁
KgCO2_beis 0 1 0.19 0.36 0 0 0 0.26 2.75 ▇▁▁▁▁

Variable type: POSIXct

skim_variable n_missing complete_rate min max median n_unique
dv_start 0 1 2022-02-11 12:00:00 2022-12-14 23:00:00 2022-07-11 05:00:00 14277

7.1.3 NG-ESO data

Figure 7.1 shows the NG-ESO half-hourly carbon intensity over time for the data period as context.

ggplot2::ggplot(elecCons_dt, aes(x = dv_date, y = dv_hms, fill = CARBON_INTENSITY)) +
  geom_tile() +
  scale_fill_continuous(name = "Carbon intensity", low = "green", high = "red") +
  labs(x = "Date",
       y = "Time of day",
       caption = "Source: NG-ESO (https://data.nationalgrideso.com/carbon-intensity1/historic-generation-mix)")
Half-hourly carbon intensity over time for the data period

Figure 7.1: Half-hourly carbon intensity over time for the data period

8 References